Stress release mechanisms for Cu on Pd(lll) in the submonolayer and monolayer regimes 
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We study the strain relaxation mechanisms of Cu on Pd(l 1 1) up to the monolayer regime using two different 
computational methodologies, basin-hopping global optimization and energy minimization with a repulsive 
bias potential. Our numerical results are consistent with experimentally observed layer-by-layer growth mode. 
However, we find that the structure of the Cu layer is not fully pseudomorphic even at low coverages. Instead, 
the Cu adsorbates forms fee and hep stacking domains, separated by partial misfit dislocations. We also estimate 
the minimum energy path and energy barriers for transitions from the ideal epitaxial state to the fcc-hcp domain 
pattern. 
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Introduction. - Metallic surfaces and nanostructures are 
essential systems for heterogeneous catalysis. Combination 
of two metals can lead to significant improvements in the 
variability and frequency of the reactions catalyzed. To pre- 
pare controlled nanostructures, it is crucial to understand the 
growth and stability of heteroepitaxial metal systems, in par- 
ticular for close-packed surfaces. On fcc(lll) it has been 
found that to release the stress, the overlayer can adopt sev- 
eral alternative strategies which lead to a structure of domains 
separated by partial misfit dislocations^"*. The domains cor- 
respond to the two favorable sites, fee and hep. This be- 
havior is expected to be ubiquitous and should not depen- 
dent strongly on the interaction potentials or the overlayer- 
substrate mismatch^Ti^. 

An interesting system, which has a large mismatch of 
—7.1%, is Cu on Pd(lll). Experimentally it has been 
found that the growth of Cu is layer-by-layer at room 
temperatureJii"— . In the submonolayer regime the substrate 
becomes gradually covered by 2D islandsii. Although near 1 
ML coverage there is some indication of dislocationsii*i^, 
only at 2-3 ML the average lattice constant reduces to the re- 
laxed value. This has been interpreted as a transition to re- 
laxed morphology ^^'^^ leading to the conclusion that the pseu- 
domorphic structure is stable up to 2-3 ML. 

In this work we unravel the atomic level stress release 
mechanisms of Cu on the Pd(l 1 1) surface from the submono- 
layer up to one monolayer regimes under conditions when 
no alloying takes place We employ two efficient computa- 
tional methodologies, namely basin-hopping global optimiza- 
tion and energy minimization with repulsive bias potential. 
Our numerical results are consistent with the experimentally 
observed layer-by-layer growth. However, we find that the 
structure of the Cu adsorbate is not fully pseudomorphic even 
at low coverages. Instead, the Cu overlayer forms fee and hep 
stacking domains, separated by partial misfit dislocations. We 
also estimate the minimum energy path and energy barriers for 
transitions from the ideal epitaxial state to the fcc-hcp domain 
pattern. 
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Figure 1: (Color online) Left panel: Low energy structure from GO 
with Ncm — 120. Pd atoms are white, while Cu colors differ accord- 
ing to Cu stacking. Right panel: The map of atomic stress in the same 
island. Stress is quantified as in Refs. [28| - [3Q| . The fee and hep areas 
are under tensile strain while the domain wall is nearly stress-free. 



Model and Methods. - In our calculations the system is de- 
scribed as a stack of rectangular (111) layers whose shorter 
edges are along close-packed rows in the x direction, longer 
edges are in the y direction and z axis is perpendicular to the 
(111) plane. Periodic boundary conditions are applied in the 
X and y directions and two bottom substrate layers are fixed to 
the bulk values. The number of atomic rows in the y direction 
is chosen to be twice that in the x direction. With this ge- 
ometry the smallest simulation cell which can accommodate 
layers with ideal lattice spacing of both Cu and Pd has 15 x 30 
edge Pd atoms. 

Metal interactions are modeled by the Embedded Atom 
Method (EAM)^^'^^, which has given a good agreement with 
experiments on the Cu/Pd(100) system^^. The energy of a Cu 
impurity in a Pd host is close to the experimental value^^ and 
the surface energies are consistent with the measurements^^- 
and first-principles studiesi^. The potential correctly predicts 
the system to be in the wetting regime, although for a bulk 
alloy it favours fee instead of the observed B2 bcc phase 
To identify the low energy configurations, we employ two 
complementary strategies, namely basin-hopping global opti- 
mization (GO) and activation-minimization (AM) techniques. 
Both methods have been described in detail elsewhere^^"^^; 
thus we will only recite their main features here. 

Global Optimization. - Basin-hopping global optimiza- 
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tion is a tool for finding, with fixed size and chemical com- 
position, the atomic configuration with the lowest potential 
energy^^.The GO tool can start the search from any given con- 
figuration, and at each step it proceeds as follows: (i) The Cu 
atoms are perturbed by the shake move described in Ref . IIO, 
and alloying is prevented, (ii) The perturbed system, including 
the Pd slab, is locally minimized, (iii) The energy of the new 
local minimum is compared to the previous one and the move 
is accepted according to the standard Metropolis criterion. 

Global optimizations for Cu island sizes A^cu > 20 start 
from an initial, randomly-arranged Cu island locally mini- 
mized on top of a 15 x 30 Pd slab with thickness varying be- 
tween 3 and 6. The two bottom layers are kept frozen in their 
ideal bulk configuration, while the other layers are allowed to 
relax during the local minimization procedure. One run of 10^ 
steps is performed for each of the systems considered. 

Activation-Minimization Procedure. - The AM technique 
starts from a perfectly pseudomorphic Cu overlayer. During 
both activation and energy minimization stages, the system 
follows the time-discretized Newtonian equations of motion 
which are solved with the standard leap-frog algorithm. We 
perform the energy minimizations with the Molecular Dynam- 
ics Cooling (MDC) technique^^. On the activation stage we 
also apply MDC but now the system is under the influence of 
the repulsive bias potential (RBP)212^. In some cases, before 
turning on the RBP, we increase the chance of getting certain 
kind of defects by perturbing some selected atoms slightly to 
a direction which is likely to initiate a nucleation process. Af- 
ter switching off the RBP we further relax the final state with 
MDC. During activation, we save intermediate states and use 
them as an input for the Nudged Elastic Band (NEB) method 
to find the minimal energy path between the initial and final 
states2&. Our final results were checked by varying the simu- 
lation box size, Pd slab thickness, perturbations and the level 
of biasing and our conclusions do not depend on these details 
either for AM or GO. 

Submonolayer regime. - In the submonolayer regime, we 
performed GO for several island sizes (5 < A^cu < 120). 
Consistently with the experimentally observed layer-by-layer 
growtbiJ^, all the global minimum energy structures for 
Ncu > 5 are flat. A low energy structure for N = 120 is 
shown in FigH] It is worth noticing that the adsorption energy 
of the island does not substantially depend on its stacking. 
When the Pd slab has stacking sequence ABC, the Cu atoms 
in the next layer can occupy sites A or B, which are referred to 
as fee or hep stackings, respectively^. The potential correctly 
assigns a lower cohesion energy to fee than to hep bulk, al- 
though the difference is small, E fee- hep < 2 meV/atom. The 
best adsorption site for a single Cu atom on Pd(l 1 1) is on the 
fee site, the energy difference from the hep stacking being 4.4 
meV. Consequently, islands with A^cu < 41 are completely on 
the fee sites. All the larger islands, where 60 < TVcu < 120, 
are in both fee and hep stackings. 

A stress color-scale graph is shown in Fig. [T] The stress of 
atom i is quantified as the isostatic local atomic pressure^"—, 
p{i) = — Tij {dE/dvij) /3 where j runs over the neigh- 
bors of atom £^ is the total potential energy of the system 
and Tij is the distance between i and j. The tensile stress is 




Figure 2: (Color online) (a) shows a stacking domain trigon in a 
hexagonal Cu island (light color) and the next-nearest substrate layer 
(dark color). When the adsorbate occupies hep sites, the dark sub- 
strate is occluded, otherwise the stacking is fcc-like and the substrate 
is visible through the overlayer. (b) shows the energy map for the 
trigon in (a). The lightly colored high energy regions are in the 
same locations as the stacking domain walls of (a) and as the highly- 
stressed region at the center of the trigon. 



largest near the center of the fee and hep domains and stress 
release originates from partial dislocations running along the 
domain walls. 

In similar systems exhibiting 2D incommensurate ordering, 
the walls are called heavy or superheavy, corresponding to 
atoms at neighbor fee and hep sites or atoms on bridge posi- 
tions, respectivel}^. The latter are more common in islands. 
In our system both walls are often distorted to bring atoms 
closer to the equilibrium distance. Especially the superheavy 
walls are typically curved and extended over several lattice 
spacings. While the wall atoms are at disadvantageous posi- 
tions, the elastic relaxation, which would not be possible with- 
out the domain structure, has longer range and releases energy 
in the whole island. Larger hexagonal islands also display the 
domains in a pattern which is called a trigon^. The center of 
the trigon shows energy increase due to high stress levels, as 
shown in FigHb). 

Monolayer regime. - Although we find that pseudomorphic 
Cu at 1 ML coverage, which is shown in FigO panel A, is sta- 
ble upon local energy minimization, both GO and AM meth- 
ods can locate lower energy structures. They consist of alter- 
nating fee and hep stacking domains similar to the submono- 
layer islands. We will refer to these structures collectively as 
the fcc-hcp domain configurations (cf. Fig. 2). We are thus 
led to investigate what is the stability of the pseudomorphic 
film relative to the fcc-hcp domains and to what extent the 
adsorbate is epitaxial. 

To mimic the row dislocation networks sometimes reported 
from similar systems^^, we perturb the 15 x 30 simulation cell 
by pushing a x directional row of the Cu layer 0.05 A down 
to the negative z direction. We activate the system with only 
3 layers of Pd, but after finding a transition path, we increase 
the substrate thickness to 9 by adding 6 layers below the orig- 
inals. As before, the two bottom of the enlarged system are 
kept frozen. The AM then yields a configuration where the 
overlayer has contracted and split to reveal an x directional 
row from the substrate, see Fig. [S] panel B. This row disloca- 
tion state has lower energy than the pseudomorphic state and 
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Figure 3: (Color online) Top: Panel A shows the 1 ML pseudomor- 
phic state, panels B and C are two fcc-hcp domain configurations 
obtained from the AM procedure. The panels A, B and C show the 
15 X 30 overlay er (light color) and the next-nearest substrate layer 
(dark color), as in Fig. [2l Bottom: Minimum energy paths connect- 
ing the pseudomorphic state to fcc-hcp domain configurations. 



contains fcc-hcp domains, but the energy barrier exceeds 1 eV. 

The fcc-hcp domains both in the submonolayer islands and 
from the GO runs at 1 ML seem to have a fairly irregular struc- 
ture. To estimate the energy barrier connecting the pseudo- 
morphic layer to these structures, we pushed 15 randomly se- 
lected Cu atoms 0.05 A to the nearest y edge direction. The 
result after the AM procedure also shows fcc-hcp domains, 
see Fig. [3l panel C. These domains have striking similarity to 
what is described in Refs. The energy path connecting 
the pseudomorphic state to the fcc-hcp domain is character- 
ized by multiple large barriers which exceed the row disloca- 
tion barrier. The holes on the overlay er are now better aligned 
to allow stress release in several directions (cf. Fig. 2, pan- 
els B and C). Because the row dislocation shows the stacking 
domains and the energy barrier is in the same range as in the 
disordered case, it is possible to view it as a special case which 
also belongs to the fcc-hcp domain class. 

The NEB energy paths of Fig. [3] suggest that if a uniform 
pseudomorphic overlay er can be formed, high energy barriers 
will prevent the decay into the fcc-hcp domain structure. If 
this overlay er was grown with conventional deposition meth- 
ods, some adatoms would eventually land on the film. There- 
fore we are led to ask if a few adatoms on a single pseudomor- 
phic layer would be enough to activate a defect nucleation. 

In our initial state we placed A/'ad Cu adatoms, all at the 
same time, on top of a pseudomorphic layer at randomly se- 
lected fee stacking positions. Substrate layer sizes 15 x 30, 
30 X 60 and 40 x 80 were used with increasing TVad- We 
find that with simple energy minimization A/'inc spontaneously 




Figure 4: (Color online) Randomly placed 0.18 ML adatoms on top 
of a pseudomorphic Cu layer incorporate spontaneously and nucle- 
ate a fcc-hcp stacking domain structure to the interfacial layer. The 
colors are as in Fig. [3 The sparse black spots are holes in the inter- 
facial layer. Insert: The average structure function in the x direction 
for eight realizations of the adsorbate layer after random incorpora- 
tion, all with the same nominal coverage 1.18 ML. The Pd-like peak, 
which is ideally above the left grey spot, is strong because the inte- 
rior of both fee and hep domains has Pd-like lattice constant. The 
peaks near the ideal Cu-position (right grey spot) are related to the 
domain walls, where most of the stress relaxation takes place. 



incorporate into the interfacial layer. For all the cases con- 
sidered A^inc saturates near the number difference between 
the fully pseudomorphic and fully relaxed layers, which is 
roughly 2/ times the layer size. The remaining Cu adatoms 
stay in the topmost layer. As a result of the incorporation, a 
fcc-hcp domain structure is again formed, as shown in FigJH 
In saturation, the fee stacking is slightly more frequent than 
hep. GO runs also display similar low-energy patterns. Over- 
layer structures with fcc-hcp domains similar to those in Fig. 
[21 have been observed on various heteroepitaxial systems on 
(111) surfaces^. 

Discussion. - The results presented here are in agreement 
with the experimentally reported layer-by-layer growth of Cu 
on Pd(lll)^^~^^, since flat islands are the low energy config- 
urations in the whole submonolayer range considered here. 
However, our results predict that the most favorable arrange- 
ment of Cu on top of Pd(l 1 1) is the fcc-hcp domain configu- 
ration rather than an epitaxial configuration. This result was 
obtained both from the global optimization method as well as 
from the activated and spontaneous relaxation starting from a 
perfect epitaxial state. 

The experimental peak intensities from Auger elec- 
tron spectroscopy (AEG), low energy electron spectroscopy 
(LEED)^^ and reflection high energy electron diffraction 
(RHEED)^^ were interpreted to indicate that the Cu deposi- 
tion is epitaxial up to 2 — 3 ML. In contrast, in our calculations 
fcc-hcp domain configurations occur already at the submono- 
layer regime. In the measured RHEED intensities, however, 
there is a Cu-like peak already visible and distinguishable 
from the Pd-like at ^ 1.5 ML^^. To estimate the degree of epi- 
taxy in the overlay ers exhibiting fcc-hcp domains, we plot the 
average structure function | ^^^^ exp(i/ca^x)p averaged over 
eight configurations after incorporation at 1.18 ML coverage 
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in FiglH The averaging is done to remove the dependence on 
the random adatom positions. Most of the signal is Pd-Hke be- 
cause both the fee and hep domains are nearly coherent with 
the substrate. However, there is some relaxation originating 
from the domain walls and this produces two peaks near the 
ideal Cu-like location, corresponding to distorted heavy and 
superheavy walls. After averaging the Cu-peaks are less vis- 
ible and for realistic temperatures the fluctuations would fur- 
ther smear the data. If the Pd signal slope were changing faster 
than the slope of the Cu signal as a function of the wave vector 
length the latter signal would be difficult to resolve. How- 
ever, there is no reason why the morphology could not change 
to completely relaxed layers from the domain structure after 
deposition of a few layers, as the experimental data shows. 
Thus, in our view, the existence of the stacking domains is not 
ruled out by the experiments. 

Similar fcc-hcp domains have been observed experimen- 
tally in Cu/Ru(0001) with misfit -5.5%i-^ Theoretical stud- 
ies suggest that the existence of the structure is not sensitive 
to details of the stressed fcc(lll) system^'^. Instead, these 
structures arise from the existence of two favorable states, fee 
and hep sites, and from the competition between strain energy 
relief and the energy cost for dislocation formation. The loss 
of epitaxy of Cu/Ru(0001) occurs above 1 ML coverage but 
in our system with misfit —7% the critical coverage is below 



1 ML, which is in line with the higher stress level. In the 
Cu/Pd(lll) system the incorporation would not lead to fur- 
ther stress release because the layer is already patterned be- 
low 1 ML. However, in the Cu/Ru(0001) system the process 
is visible. The XPS intensity data for Cu/Pd(l 1 1) system also 
leave room for layer-by-layer morphology above 1 MLi^. 

Conclusion. - We have used computational approach to un- 
derstand the structure of Cu films on the Pd(lll) surface in 
the low coverage regime extending from 20-atom islands to 
1 ML. Our results agree with layer-by-layer growth, as is ob- 
served experimentally. Our model predicts that the epitaxi- 
ality is lost already below 1 ML coverage and Cu forms hep 
and fee stacking domains, separated by wall regions in which 
Cu recovers its lattice spacing and releases its tensile stress. 
This pattern, similar to what is seen in Cu/Ru(0001) at higher 
coverages, encourages further experimental investigations di- 
rectly probing the overlay er structure. 
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